%%% Written by  Daniel Lewis, 2020, for "Robust Inference in Models 
%%% Identified via Heteroskedasticity", Review of Economics and Statistics.
%
% For a 2-variable model, this function takes as inputs the Tr by 3
% matrices of the unique vectorization of outer products of residuals for 
% each regime, the total sample length, and the parameter vector. It 
% returns the GMM objective function using the efficient weighting matrix,
% the sum of moments, the covariance of moments, and the efficient 
% weighting matrix.
%
% Inputs:
% etasq1:   outer product of regime 1 residuals (T1x3 unique vectorization)
% etasq2:   outer product of regime 2 residuals (T2x3 unique vectorization)
% T:        total sample length
% p:        parameter vector
% Outputs:
% obj:      GMM objective function
% f_T:      sum of moments 
% omega:    moment covariance matrix
% matinv:   efficient weighting matrix

function [obj,f_T,omega,matinv]=momfast(etasq1,etasq2,T,p)
V=[1,p(2);p(1),1]*diag([p(3),p(4)])*[1,p(2);p(1),1]'; % construct fitted reduced form covariance matrix for regime 1
V1=[V(1:2),V(4)]; % vectorize unique elements
V=[1,p(2);p(1),1]*diag([p(5),p(6)])*[1,p(2);p(1),1]'; % construct fitted reduced form covariance matrix for regime 2
V2=[V(1:2),V(4)]; % vectorize unique elements
f_T=[[etasq1-V1;zeros(length(etasq2),3)] [zeros(length(etasq1),3); etasq2-V2]]; % compute moments and stack
omega=f_T'*f_T/T; % compute covariance matrix of moment equations
f_T=sum(f_T)'; % compute sum of moments
matinv=nearestSPD(omega^-1); % compute efficient weighting matrix (nearestSPD corrects in case of numerical issues with symmetry or positive definiteness)
obj=f_T'*matinv*f_T/T; % compute GMM objective function
end
